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We present a non-canonically symplectic integration scheme tailored to numerically computing 
the post-Newtonian motion of a spinning black-hole binary. Using a splitting approach we combine 
the flows of orbital and spin contributions. In the context of the splitting, it is possible to integrate 
the individual terms of the spin-orbit and spin-spin Hamiltonians analytically, exploiting the special 
structure of the underlying equations of motion. The outcome is a symplectic, time-reversible 
integrator, which can be raised to arbitrary order by composition. A fourth-order version is shown to 
give excellent behavior concerning error growth and conservation of energy and angular momentum 
in long-term simulations. Favorable properties of the integrator are retained in the presence of weak 
dissipative forces due to radiation damping in the full post-Newtonian equations. 
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I. INTRODUCTION 

In general relativity there is a rich phenomenology as- 
sociated with two orbiting compact objects and the grav- 
itational waves that are emitted in the process. These 
phenomena can be studied to very high accuracy in the 
post-Newtonian approximation of general relativity, see 
e.g. pp. The post-Newtonian equations of motion are 
well-known up to a certain order, and numerical integra- 
tion of the resulting ordinary differential equations can 
be performed to obtain the orbits without immediate dif- 
ficulty, at least as long as certain evolution times and 
accuracy requirements are not exceeded. 

In this paper we address the question of long-term in- 
tegration of the post- Newtonian equations of motion, al- 
lowing unequal masses and spins. A study of a large num- 
ber of orbits is of interest for various reasons, of which 
we want to mention only one. In general relativity, there 
is the possibility of chaotic orbits, which can place severe 
demands on the quality of the numerical integrator. 

The occurrence of chaotic trajectories of orbiting spin- 
ning binaries has been discussed, for example, in [2HS]j 
implying unpredictable irregularities in the waveforms. 
However, there still seems to be a controversy among dif- 
ferent authors under what conditions the post-Newtonian 
equations lead to chaotic motion. Various possible indi- 
cators for chaos like Lyapunov exponents, Poincare sec- 
tions, and basin boundary plots, are employed in the in- 
vestigations found in the literature, and interesting re- 
gions of the configuration space, including mass ratio, 
eccentricity, spin orientations and lengths, are densely 
covered by thousands of numerical simulations. Typical 
studies of chaos require numerical long-term evolutions of 
the system, depending on the timescale at which chaotic 
behavior becomes apparent. 

The structure of the post-Newtonian equations of mo- 
tion permit a separate investigation of bodies orbiting 
in a conservative fashion (not incorporating radiation 
damping) and of inspiraling masses losing energy by emis- 
sion of gravitational waves. The conservative case fa- 



cilitates the detection of chaotic tendencies as simula- 
tions can in principle be carried out arbitrarily long, 
whereas in the dissipative scenario the objects would 
finally merge. Both kinds of simulations, conservative 
and non-conservative, would benefit from an efficient and 
well-behaved long-term integrator. 

Computational challenges from a wide variety of areas 
of research like celestial mechanics, molecular dynamics, 
quantum mechanics, as well as abstract numerical anal- 
ysis motivated the development of structure-preserving 
algorithms for differential equations. Geometric integra- 
tors such as symplectic or symmetric methods are known, 
from numerical experience and theory, to yield a substan- 
tially improved long-term behavior compared to standard 
integrators such as explicit Runge-Kutta methods. The 
improved behavior concerns the exact or approximate 
preservation of conserved quantities (like energy and an- 
gular momentum) without drift, slower error growth, and 
also a more faithful representation of Poincare sections 
as used in investigations of chaotic behavior. We refer to 
the monographs [7J [5] and numerous references therein. 

In the following we construct a structure-preserving, 
efficiently implementable integrator for the equations 
of motion of binary spinning black holes in the post- 
Newtonian approximation. The integrator is based on 
a splitting into conservative motion and dissipative per- 
turbation due to radiation-reaction forces, into orbital 
and spin evolution, and on further splittings between 
the Newtonian and post-Newtonian Hamiltonian and be- 
tween different terms in the spin Hamiltonian. 

In the conservative case, the resulting algorithm is a 
Poisson integrator (or non-canonically symplectic inte- 
grator in another terminology) and time-reversible. The 
term 'non-canonical' refers to the spin algebra, which 
yields a non-canonical Poisson bracket for which no stan- 
dard structure-preserving integrators exist. The Poisson 
integrator preserves the spin lengths exactly. Our nu- 
merical experiments show no drift in energy and total 
angular momentum. 
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The paper is organized as follows. In Section II 
we present the post-Newtonian Hamiltonian equations 
of motion in ADMTT gauge (Arnowitt-Deser-Misner 
transverse-traceless) . The construction of the symplectic 
splitting integrator, which exploits the special structure 
of that Hamiltonian, is described in Section III. For the 
conservative equations the performance of the new in- 
tegrator is contrasted with the behavior of the classical 
Runge-Kutta method in Section IV. Finally, in Section V 
we incorporate the dissipative terms into the integrator. 



with order. However, the complexity of the single terms, 
whose explicit form can be found in [9Hl2j , increases con- 
siderably with order, which seriously affects the compu- 
tational costs. 

The spin contribution in Eq. Q also forms a series 
expansion. We are going to consider the leading-order 
terms proportional to G/c 2 and up to quadratic in the 
spins only. These expressions can be grouped in spin- 
orbit interactions and spin-spin interactions |131 114] , 



H 



Spin 



(6) 



II. POST-NEWTONIAN EQUATIONS OF 
MOTION FOR SPINNING BINARY SYSTEMS 

This section presents the equations of motion of 
a black-hole binary system consisting of objects with 
masses m a , positions X a , momenta P a and spins S a 
(a = 1,2). For our purposes it is sufficient to restrict 
considerations to the center-of-mass dynamics, where 

P = Pi= P 2 

The Hamiltonian formalism turns out to be very use- 
ful when working with the canonically conjugate position 
and momentum variables and in distinguishing conser- 
vative and radiative effects. The system's equations of 
motion take the form 



dX -fx m- dH 



dp 

lit 

dS a 



= {P,H}+F = - 



dH 
dX 



F, 



(i) 

(2) 
(3) 



where F is the non-conservative force and x denotes the 
usual vector cross product. The conservative part con- 
sists of an orbital and a spin contribution with the com- 
posed Hamiltonian 

H(X, P, S 1 ,S 2 ) =H OTh (X, P) + H Spin {X, P, S U S 2 ). 

(4) 

Depending on the choice of gauge, different formulations 
are possible. We use the 3PN accurate orbital Hamilto- 
nian in ADMTT gauge derived by Damour, Jaranowski, 
Schafer [5HT2"]. It forms an expansion in the parameter 
1/c 2 and goes beyond the classical Newtonian Hamilto- 
nian £/n by adding post-Newtonian, relativistic correc- 
tions, 

H OTh (X, P) = Mc 2 + H N (X, P) + \h 1pn {X, P) 

Cr 

+ 4#2PN(X, P) + 4#3PN(X, P), 

(5) 

where M = mi + m 2 denotes the total rest mass. Al- 
though the convergence of this series can be slow, the 
magnitude of the energy and force contributions asso- 
ciated with the individual terms nevertheless decreases 
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mi ' 
(10) 

(11) 



Here, L = X x P is the orbital angular momentum and 
N is the unit vector X/R, where R= \X\. We assume 
the Newton- Wigner spin supplementary condition, which 
provides the notion of the spin vector and affects the 
definition of the bodies' worldlines. 

Turning towards the radiation-reaction-force term F 
to be applied in Eq. ([2|, we use the expression derived 
by |15) in cases where we leave the conservative regime. 
The expression reads 



1 dE 
"uAL\ ~dt J 



V n 

15 L 2 R 



Gl 



48^V-Si 
mi J 



61 + 48 



mi 
m 2 



P-S 2 >L, 



(12) 



where v — fi/M with fj, = m\m 2 jM being the reduced 
mass. It is crucial to note that the force expression de- 
pends on the orbital angular frequency ui, which is also 
hidden in the invariant velocity parameter v u , 



/GMw\ 1/3 



(13) 



The problem with u) is that it brings the velocity X = 
dX/dt into play, 



\X - N(N ■ X)\ 
R ' 



(14) 



which is inconvenient when dealing with momenta. Note 
that the link between the momentum P and the veloc- 
ity is given by the complicated post-Newtonian relation 
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X = dH{X,P,Si,S 2 )/dP. The parameter lj is also 
needed to express the energy-loss rate dE/dt appearing 
in Eq. (Tl2|). We will use the 3.5PN accurate flux for circu- 



larly orbiting masses which has been derived by 
and is given in [TS] as 



dE 
~dt 



32 



i + h{y)< + [hi") + hsoK 



+ [Uiy) + UssH + h(»)v 5 u + hM 
+ f m v% hx{Av u ) + f 7 {v)vl 



(15) 



The coefficients of this series are constant, except for / 3 so 
and /4ssj which depend on the spins and the direction of 
L. 

It is worth stressing that the energy flux in this form is 
strictly valid only for planar, circular motion without any 
trace of eccentricity. The spin-orbit interaction, Hso, is 
known to produce a precession of the orbital plane. While 
this slow precession is supposed to give only higher-order 
deviations in the energy flux, the orbital oscillations pro- 
voked by the spin-spin terms would lead to noticeable in- 
accuracies. A more complicated version of the flux which 
is valid also for eccentric orbits has been derived by [TS] • 

However, for our purposes the dissipation according to 
Eq. (151 is considered sufficiently accurate. We do not 



consider it as the goal of this work to evaluate the range of 
validity of the equations, such as the restriction of Eq. Q 
to small spin magnitudes, for example, or the complete 
breakdown of the equations of motions at small binary 
separations. Instead we want to take the equations, or 
parts of them, as given and present numerical methods 
well-suited for these kinds of equations. 



III. A SYMPLECTIC SPLITTING 
INTEGRATOR FOR THE CONSERVATIVE 
SYSTEM 

We first focus on the conservative motion, with the 
dissipative force F set to zero. The system 0-(|3]) is then 
a non-canonical Hamiltonian system (or Poisson system) 
with the Poisson bracket composed of a canonical bracket 
and the two spin brackets, 



{F,G} 



/ dF dG 
\d~X ' dP 

2 

+ ^det 



dF dG 
~ dP ' dX 

dF dG 



dS n 



8S n 



(16) 



The exact flow : (X(0), P(0), Si(0), S 2 (0)) h> 

{X(t), P(<), Si(t), S2{t)) is a Poisson map (or non- 
canonically symplectic map in another terminology) : For 
all smooth functions F — F(X , P, Si, S2) and G, the 
flow preserves the bracket as 



{F olf ?,Go^} = {F,G}^ 



ft 



H 



(17) 



In addition to the total energy H, the spin lengths \S a \ 
and the total angular momentum J = L + Si + S2 are 
conserved quantities. 

For the numerical treatment of this system we pro- 
pose a splitting integrator with the following structure- 
preserving properties (see, e.g., [Jj for terminology): 

• The method is a Poisson integrator (or non- 
canonically symplectic integrator) for the bracket 
( |16[ ): For all smooth functions F and G, the 
discrete flow over a stepsize h preserves the 
bracket, 



{Fo$«Go$f} = {F,G}o$f. 



(18) 



The method preserves the spin lengths \S a \ (the 
Casimirs of the bracket). These properties imply 
that a step of the method equals the exact flow 
(up to terms that are exponentially small in the in- 
verse stepsize) of a Poisson system with the original 
bracket ( |16[ ) and a slightly modified Hamiltonian 
(see [7J Chap. IX]). 

• The method is symmetric and therefore preserves 
all reversal symmetries present in the system. 

The basic method is of second order, but it can be en- 
hanced to higher orders by composition. 

The method is based on a splitting of the Hamiltonian 
into its orbital, spin-orbit and spin-spin contributions, 



H = H, 



Orb 



SO 



H SSl 



(19) 



with further splittings of -Horb, Hso an d -ffss to be de- 
scribed below. Accordingly, we start from the approxi- 
mation of the flow of the system |lj"([3]) with F = 
over a time step h by the symmetric splitting 



(20) 



The individual flows in this formula will now be further 
approximated in a structure-preserving way. 

Orbital integrator. For the approximation of the 
flow of the orbital Hamiltonian -fforb = -^n + -HpN we 
further split into the Newtonian part, which yields pure 
Kepler motion, and the computationally expensive Post- 
Newtonian corrections: 



#PN 



0( Ph/2- 



(21) 



Besides the possibility to use the exact Kepler flow tp^ 
we employ a high-order symplectic approximation to 
it. The sixth-order method labeled p9s9 obtained by 
symmetric composition of Stormer-Verlet steps (see 
Sect. V. 3. 2]) showed best performance in our case. 

The post-Newtonian flow </?^/ p 2 N is approximated by the 
symplectic Euler scheme 



h/2 ■ ra+1/2 A n + ^ p, D n+1/2 , f 



X n+ i/2 -X n , 2 dp 

P n +l/2 =P n - ^ { X n+l/2,Pn) (22) 
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and its adjoint (exchange non+l and h <H> —h), 

*/i/2 • - A «+l/2 + 2 Qp ( A «+l/2,-^«+lJ 

p p ?i d g PN 

*7»+l — r n+1/2 — 2 QX \ n+1 / 21 n+1 ) 



(23) 



For the non-separable post-Newtonian Hamiltonian i/pN, 
these steps are implicit in X n+1 / 2 for 5? h and in P n +i 

for QfJ^'* ■ The post-Newtonian corrections are expected 
to be only minor, and so the initial values of -X" and P 
are excellent starting values for the fixed-point iterations. 
Our orbital integrator thus approximates the flow <^° rb 
by 



h/2 



o $f N o $ 



h/2 



Rotations. For a Hamiltonian 



H r 



with a constant vector 17, the equations of motion 

S a = J~i X S'a 



(24) 



(25) 



(26) 



are readily solved by a rotation given by Rodrigues' for- 
mula, 



s a (t) = n(n,t) s B (o) 



(27) 



S„(0)+«^x Sa (0) 
1 /sin(^|f2|)\ 2 



Alternatively, the rotation can be efficiently implemented 
using quaternions (see, e.g., Sect. VII. 5. 3]). No con- 
stant rotation vector f2 appears in the spin-orbit and 
spin-spin Hamiltonians ([7])- (10). However, applying fur- 
ther splittings and appropriate reformulations will nev- 
ertheless allow us to make use of the rotation formula 
for each of them. We note in passing that although for 
the full system the spin evolution cannot be determined 
analytically, a perturbative calculation exists for approx- 
imately equal masses for the spin-orbit part [20 . 

Spin-orbit integrator. Up to a constant factor, 
which we omit in the following, the spin-orbit Hamil- 
tonian ([7]) is given by 



Hso(X, P, Si, S 2 ) — S c tf ■ L/R 3 



(28) 



with S e g = C1S1 + c 2 S 2 . A way to arrive at a symplectic 
approximation to its flow, is to split 

H so = Hlo + Hlo + Hlo with Hi Q = Sl s V/R 3 . (29) 

As we show next, the equations of motion for Hg can 
be solved exactly by rotations of X, P, Si, S 2 . Denoting 



by e, the ith unit vector, the equations of motion read 
(for i = 1, and analogously for i = 2,3) 

X = Sl sei x X/R 3 (30) 
P = Slg ei x P/R 3 + 3^o X/R 2 (31) 
S a = c a L 1 e 1 x S a /R 3 . (32) 

The first equation implies ^R 2 = 2X • X = 0, so that 
R = const. The spin equations imply = and hence 
S^ff = const. This shows that the equation for X is solved 
by the rotation lZ(Sl s ei/R 3 , t). The equation for P is 
then also solved analytically. For the angular momentum 
we obtain from the equations of motion for X and P the 
differential equation 



L = S 1 cS e 1 xL/R 3 , 



(33) 



which implies L 1 = 0, so that also L 1 — const., and hence 
the spin equations are also solved by simple rotations. 

H 1 

This way, the flow <p t so is obtained by four rotations: 

X(t) = K(Sl s e 1 /R 3 ,t)X(0) (34) 
P(t) = K{Sl sei /R 3 ,t) (P(0) + ZtHl o X{0)/R 2 ) (35) 
S a (t)=K(c a L 1 e l /R 3 ,t)S a (0). (36) 

We then approximate the flow ^J 2 by the splitting 



4> 



Hi 



Hi 



Hi 



h/2 = <P h /2 ° Vh/2 ° <Ph/2 



(37) 



or by the adjoint method, with reversed order of the 
flows, 



Hi 



Hi 



Hi 



h/2 = <PhJ* ° ° 



(38) 



Spin-spin integrator. In the following an approx- 
imation to the spin-spin coupling flow ip h f 2 is derived 
by a splitting of the spin-spin Hamiltonian, where again 
each part is readily integrated exactly by rotations. For 
related splitting approaches for classical spin systems we 
refer to [2lti23] . see also [24]. 

Each of the terms appearing in the spin-spin Hamilto- 
nian (16| belongs to one of the forms (omitting constant 
factors) 



H A = S x ■ S 2 /R 3 , 

H B = S a ■ S a /R 3 - 

H c = (S 1 -N)(S 2 -N)/R 3 , 

H D = h 2 (S a -N) 2 /R 3 , 



(39a) 
(39b) 
(39c) 
(39d) 



Since there is no dependence on P in any of these Hamil- 
tonians, it follows that X(t) remains constant: 



X(t)=X(0). 



(40) 



The dependence on X is through the factor 1/R 3 and 
in cases C and D also via the unit vector N = X/R. 
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In case A we have dH A /dX = -3H A X/R 2 , which is 
constant along the evolution. Therefore, we obtain 



P(t) = P(Q) + MH A X/R 2 



(41) 



The analogous formula is valid also for H B . Due to the 
cross product in ((3j|, the Hamiltonian H B does not con- 
tribute at all to the motion of the spins. This feature can 
be exploited for the integration of the spin equations for 
H A , which read 

St = S 2 x St/R 3 = {St + S 2 ) x St/R 3 , (42) 
S 2 = St x S 2 /R 3 = (St + S a ) x S 2 /R 3 , (43) 

where consequently the sum St + S 2 is found to be zero. 



Hence, St + S 2 must be constant and Eqs. (42 1 and (43) 
are therefore solved by rotation, 



s a (t) = n((St(o) + s 2 (o))/R 3 ,t) s a (o). 

The Hamiltonian H has the spin equations 

St = (S 2 - N)N x St/R 3 , 
S 2 = (St-N)Nx S 2 /R 3 . 



(44) 

(45) 
(46) 



Taking the inner product with TV shows that both St ■ N 
and S 2 • N are constant. The equations are therefore 
solved by rotation, 

St(t) = K((S 2 (0) ■ N)N/R 3 ,t) St(0) (47) 
S 2 (t) = U((St(0) ■ N)N/R 3 ,t) S 2 (0). (48) 

For H , the equations of motion for the momenta read 

P = 5H C X/R 2 -S 1 (S 2 -N)/R i -S 2 (SfN)/R' i , (49) 

which are solved by direct integration: 

P(t) =P(0) + 5tH c X/R 2 (50) 

- Q((S 2 (0) ■ N)N/R 3 ,t)St(0)(S 2 (0) ■ N)/R i 

- Q((St(0) ■ N)N/R 3 ,t)S 2 (0)(St(0) ■ N)/R 4 , 

where 

C(n,t)s a (o)= / n(n,t')s a (o)dt' (51) 

Jo 



tS a (0) 



1 /sin(|t|n| 



n x s a (o) 



t\fl\ - sin(i|fi|) 

W 3 



n x n x s a (o). 



The treatment of the Hamiltonian H D is similar to the 
just expounded procedure and is therefore not shown. 
Taking all parts together, we approximate the spin-spin 
flow ^J 2 by the composition of the exact flows of the 
subsystems, written schematically as 



Its adjoint ® h J 2 * features the reversed order of all con- 
tributions. We have chosen an ordering of the vari- 
ous sub-Hamiltonians which might not be optimal. We 
note, however, that the flows of A and B commute, since 
{H A ,H B } = 0. This implies that no splitting error oc- 
curs between A and B. Similarly, the flows of C and D 
also commute. As compositions of exact flows, the dis- 
crete spin-spin flows and ^f 2 '* are Poisson inte- 
grators. 

Higher order by composition. The overall method 



4> 



h/2 



o $ 



^SO 1 

h/2 



o $f °' b o $f ?° $Hss (53) 



h/2 



h/2 



is a Poisson integrator as a composition of Poisson maps. 
It is only of second order, but its order can be enhanced 
by composition; see, e.g., [7] Chap. V. 3] or In our 

numerical experiments we will lift its order to four by 
Suzuki's composition scheme |26| : 



H , 4th order 



$ H \ o $ ff . o $ ff , o $ H h o $ H , 

■ysh 74 n 73" 72^ 71 ft' 



where the lengths of the substeps are 

1 



41/3 



7i = 72 = 74 = 75 = 



4 1 / 3 



73 = 



41/3' 



(54) 



(55) 



The just introduced integration scheme Eq. (53) shall 



be termed SPN integrator (symplectic post-Newtonian 
integrator). The fourth-order version Eq. (54) will be 
referred to as SPN4 in the following. 



IV. NUMERICAL EXPERIMENTS 

For the demonstration of the properties and perfor- 
mance of the symplectic integrator we use the follow- 
ing configuration, consisting of two black holes with non- 
trivial, maximal spin vectors: 



mt = 0.25 
m 2 = 0.75 

R = (50,0,0) 

P = (0, 0.027475637, 0) 



(56) 



St = 

s 2 = 



(-1,0,0) 

(1/V2, 0, 



1 



These initial data yield low-eccentricity orbits (see |27j). 
the average period of which is approximately T« 2, 286. 

For comparison with the newly developed fourth-order 
symplectic integrator we use the classical fourth-order 
Runge-Kutta method (see, e.g., [28]) in a self-made C 
version. Inadequacies in the efficiency of the implemen- 
tation of either the RK method or the symplectic scheme 
would only affect the CPU timings of the subsequent re- 
sults but not the principal behavior and insights. 
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FIG. 1: Comparison of error growth between SPN4 and RK4 
method for different step sizes. 



We measure the accuracy of the numerical integration 
by the scaled error norm 



fAXY /APY 



ASi 



AS 2 



1/2 



(57) 



where the individual terms are computed as the Eu- 
clidean norm of the componentwise deviation from a ref- 
erence solution (RK4, step size h = 0.0125), for instance, 



AX 

IT 



E 



Xj- X. 
Xf 



ref 



1/2 



(58) 



In the SPN4 method there is a fixed-point iteration 



to solve the implicit equations (22|-(23|. As a stopping 



criterion for this iteration we use the condition 

1/2 



3 

E 

,i=i 



XV 



xr 1 ) 2 1 



3 

E 

i=l 



{xi 



for Eq. (22), and 



' 3 

E 

.i=l 



PT 



pr 1 ) 2 1 E(o 2 
i=i 



1/2 



< £fpi (59) 



< cfpi (60) 



for Eq. ( 23 1 . The number of iterations will be limited 



to iVppi. In the comparisons below we chose the com- 
bination eppi = I0~ 12 and iVppi = 5, which showed ap- 
propriate performance for our precision demands. The 
step size h is chosen to give roughly comparable com- 
puter time for the two methods, implying a larger h for 
SPN4 in our examples. 



As indicated in Fig. [T] the error growth of the symmet- 
ric, symplectic method SPN4 is linear, a feature which is 
not unexpected (cf. [7J Chaps. X and XI]). By contrast, 
the non-symmetric, non-symplectic RK method suffers 
from quadratic error growth, making it inferior for long- 
term simulations. In the example the simulation time is 
t = 512,000 which gives about 224 revolutions of the bi- 
nary black-hole system. Supposing approximately equal 
numerical costs, there evidently exists a time ibrcak even 
beyond which the SPN4 method gives more accurate so- 
lutions than RK4. An analysis of the CPU time spent for 
the different Hamiltonians in the SPN4 method reveals 
that about 88% is consumed by the ffpN part, 9% by 
i?N, an d only 3% for all the spin parts together. 

Fig. [T] also depicts the 4th-order convergence of both 
methods, SPN4 and RK4. Dividing the step size by two 
yields an error which is smaller by a factor 1/16. 




0.2 0.3 
t I 10 6 

FIG. 2: The symplectic method (h = 64) preserves a nearby 
Hamiltonian, whereas the RK integrator (h = 8) produces a 
drift. The system's exact energy is about Eq w —1.817- 10 -3 . 




0.01 0.02 0.05 0.1 0.2 
t I 10 6 

FIG. 3: The variations of J induced by the symplectic method 
can be attributed to machine-precision effects. 

Another striking feature of the symplectic integrator is 
the near-conservation of the energy of the system. Fig. [2] 
shows that only minor oscillations in the system's en- 
ergy occur, whereas a linear drift results for the non- 



7 



symplectic RK method, becoming severe for long-term 
simulations. A similar observation is made for the other 
conserved quantities J and S a , that is the absolutes of 
the total angular momentum and the spins, cf. Figs. [3] 
and|U 




t I W 



FIG. 4: The spin lengths are preserved exactly by the SPN 
method. Oscillations are due to machine-precision effects. 
The RK integrator yields a drift. 



V. ADDING DISSIPATION 

So far the non-conservative term F in Eq. ^ has 
been omitted, and the construction of the SPN method 
was targeted to the conservative equations. However, in 
the physically realistic scenario dissipation has to be ac- 
counted for. Energy is transferred into propagating cur- 
vature of spacetime, the radiation of gravitational waves. 

As is known from Chap. XII], it is advantageous 
for dissipatively perturbed systems to use numerical 
methods that reduce to a symmetric and/or symplec- 
tic method for the unperturbed, conservative system. 
Therefore we use the SPN method as derived so far and 
incorporate the dissipative force F via a further splitting. 
The starting point is the symplectic flow before it is lifted 
from second to fourth order by the Suzuki composition 
scheme (54 1. The dissipative contribution, denoted by 



i fih/2 , > * s P a dded around the conservative flow, 



,„diss ,F „ ,„cons „ ,F 
Vh ~ Vh/2 ° Vh ° <Ph/2- 



(61) 



In the numerical implementation the conservative flow 
is again approximated by Eq. ( 53 1 , and the dissipative 



perturbation is modeled by an Euler step 



3 



-F(X 7 P } X) with X = ^ p (X,P) 



(62) 



The adjoint method is the implicit Euler method, 



$* h F : P 



ll ~ ' ~ pjTT _ 

-F(X,P,X) with X = _(X,P). 

(63) 



This results in a symmetric, second-order integrator 
given by 



jjjCUSS _ 



h/2 



Note that we have to provide the computationally ex- 
pensive term dH/dP in order to compute the force F. 
Retaining the symmetry of Eq. (64) is, however, of no 



structural importance for the dissipative contributions. 
Hence, if we replace the implicit Euler method by the 
explicit method resulting from two fixed-point iterations, 
we obtain overall a second-order method that uses three 
additional evaluations of dH/dP to incorporate the dis- 
sipative forces. Applying Suzuki composition, Eq. (54 1, 
the order of the method then is increased to four 
higher) . 



oi- 




lO- 16 

0.01 0.02 



0.05 



0.1 



0.2 



0.5 



t/t n 



Method -Ratart h *cpu 

SPN4 200 256 6.1s 

SPN4 100 128 6.3 s 

SPN4 50 64 6.6 s 
Linear growth 



Method fl H tart h t C pu 
RK4 200 32 10.3 s 
RK4 100 16 10.1s 
P..K4 50 8 10.3 s 
Quadratic growth 



FIG. 5: Error growth and energy deviation for the dissipative 
implementation of SPN4 and the RK4 method starting from 
different initial separations i? s tart- Applying different step 
sizes while keeping the number of steps (8000 for SPN4, 64000 



for RK4) results in different simulation times t n 
normalization of the time axis. 



used for 



For the numerical experiments the initial data Eq. ( 56 ) 



are modified to still give low-eccentricity orbits in the 
dissipative regime. For example, an inspiral starting at 
R = 50M would have the initial momentum 



P = (-3.5267394 • KT 6 , 0.027475637, 0), 



(65) 



obtained with the algorithm presented in [27]. The en- 
ergy loss in the system now causes the separation of the 
black holes to decline from R — 50 down to R = 36.7 
during the simulation. In general, there would be a need 
for stepsize control mechanisms, especially for longterm 
simulations. A reversible step size strategy as in [71 
Sect. VIII. 3. 2] can be employed based on the conservative 
part, but this is beyond the scope of the present work. In 
our numerical experiments, a constant stepsize is chosen 
such that sufficient resolution is granted throughout the 
simulation. 

Variation of the initial separation i? s tart provides a 
possibility to adjust the degree of dissipation. For a bi- 
nary separation of 50M the perturbation is already quite 
strong. The ratio of dissipative force to conservative 
terms amounts to \F\/\dH/dX\ « 6.0- 10~ 5 . For separa- 
tions of 100M and 200A/ the ratio goes down to 1.1 ■ 10~ 5 
and 2.0 • 10 -6 , respectively. Fig. [5] shows how increas- 
ing dissipation affects the performance of the integrators. 
The linear error growth of the SPN4 integrator does not 
persist for strong dissipation. However, its performance 
remains superior to the RK4 method, which also detaches 
from the quadratic error growth. In a similar fashion the 
ideal behavior of the SPN4 method concerning energy 
and angular momentum, i.e. minor fluctuations around 
the exact values, is lost. Even for very small pertur- 
bation with i? s tart = 200M we observe a linear growth 
exceeding the machine noise. For even stronger dissipa- 
tion the growth turns out to be faster than linear. The 
very same behavior is seen for the angular momentum 



(not shown). The spin lengths still remain preserved by 
the SPN4 method, due to the construction of the algo- 
rithm. In our numerical experiments the SPN4 method 
still shows a better performance than the RK4 method 
for sufficiently long integration intervals. 



VI. CONCLUSION 

We presented a new numerical scheme designed for in- 
tegrating the post-Newtonian equations of motion for a 
spinning black-hole binary. The special structure of the 
Hamiltonian enabled us to use a splitting algorithm that 
exactly preserves the non-canonical symplectic structure 
of the system. The result is an integrator with excellent 
long-term performance. The integrator was also modi- 
fied to include dissipative effects contained in the more 
realistic 3.5PN accurate model. Even though symplec- 
ticity is no longer present in the dissipative case, the 
integrator was shown to exhibit favorable properties in 
long-term simulations with small dissipative forces. With 
these properties we expect this new integrator to be a 
useful and efficient tool in the investigation of chaotic 
traits in the dynamics of spinning binary systems, signif- 
icantly reducing the risk of numerical artifacts compared 
to standard integrators. 
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